main.R

#:::::::::::::::::::::::NORMAL DISTRIBUTION::::::::::::::::::::::::::::::::::::#
# ******** THE PROBABILITY DENSITY FUNCTION ********:

x= seq(-4, 4, length=200)
y= 1/sqrt(2*pi)*exp(-x^2/2)
plot(x, y, type="l", lwd=2, col="red") 

# is the same as:

x = seq(-4, 4, length=200);
y = dnorm(x, mean=0, sd=1);
plot(x, y, type="l", lwd=2, col="red");

# Q: Check the effect of changing the mean and then the standard 
# deviation and notice the effects onthe normal distribution plots. 
# What kind of parameters are these?

# Increasing mean:
for (i in 1:5)
{
  x = seq(-4, 4, length=200);
  y = dnorm(x, mean=i, sd=1);
  plot(x, y, type="l", lwd=2, col="red");
}

# Increasing Standard Deviation:
i = 0;
for (i in 1:5)
{
  x = seq(-4, 4, length=200);
  y = dnorm(x, mean=0, sd=i);
  plot(x, y, type="l", lwd=2, col="red");
}

# Increasing Mean & Standard Deviation:
i = 0;
for (i in 1:5)
{
  x = seq(-4, 4, length=200);
  y = dnorm(x, mean=i, sd=i);
  plot(x, y, type="l", lwd=2, col="red");
}

# A: It is obvious that changing the mean will result in a shift of the graph,
# and changing the standard deviation will result in a scale of the graph.
# Changing both the mean and standard deviation will result in a shift and scale.
# The mean is the location parameter, and the variance is the shape parameter.

# ******** 68%-95%-99.7% Rule {Empirical Rule} ********:

# 1 SD From mean:
x=seq(-4,4,length=200)
y=dnorm(x)
plot(x,y,type="l", lwd=2, col="blue")
x=seq(-1,1,length=100)
y=dnorm(x)
polygon(c(-1,x,1),c(0,y,0),col="gray") 

pnorm(1,mean=0,sd=1)-pnorm(-1,mean=0,sd=1) 

## [1] 0.6826895

# 2 SD From mean:
x=seq(-4,4,length=200)
y=dnorm(x)
plot(x,y,type="l", lwd=2, col="blue")
x=seq(-2,2,length=100)
y=dnorm(x)
polygon(c(-2,x,2),c(0,y,0),col="gray") 

pnorm(2,mean=0,sd=1)-pnorm(-2,mean=0,sd=1) 

## [1] 0.9544997

# 3 SD From mean:
x=seq(-4,4,length=200)
y=dnorm(x)
plot(x,y,type="l", lwd=2, col="blue")
x=seq(-3,3,length=100)
y=dnorm(x)
polygon(c(-3,x,3),c(0,y,0),col="gray") 

pnorm(3,mean=0,sd=1)-pnorm(-3,mean=0,sd=1) 

## [1] 0.9973002

# ******** TESTING QUANTILES ********:
qnorm(0.95, mean=0, sd=1)

## [1] 1.644854

#Hence, there is a 95% probability that a random number less than or equal to
# 1.645 is chosen from the standard normal distribution. 

qnorm(0.20, mean=0, sd=1)

## [1] -0.8416212

# We now know that the probability of selecting a number from the standard 
# normal distribution that is greater than or equal to -0.842 is 0.80. 

# ******** DRAWING A RANDOM SAMPLE ********:
sample1 = rnorm(5000,10,4)
summary(sample1)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  -3.998   7.214  10.010   9.965  12.710  24.530

var(sample1)

## [1] 16.08548

# Q: Do these statistics match the parameters?
# A: Yes. We see that the mean is around 9.911, which is close to 10. We also see
# that sigma^2 = variance is around 16.16, which is close to sigma = 4.

#:::::::::::::::::::::::END NORMAL DISTRIBUTION:::::::::::::::::::::::::::::::::#

# PART (A):
#:::::::::::::::::::::::1-GAMMA DISTRIBUTION::::::::::::::::::::::::::::::::::::#
# ******** THE PROBABILITY DENSITY FUNCTION ********:

x= seq(0, 10, length=200)
alpha = 5; 
lamda = 2;
y= (x^(alpha-1)*exp(-lamda*x)*(lamda)^alpha / gamma(alpha));
plot(x, y, type="l", lwd=2, col="red") 

# is the same as:

x = seq(0, 10, length=200);
y=dgamma(x,5,2) 
plot(x, y, type="l", lwd=2, col="red");

# Q: Check the effect of changing the two parameters and notice the effects on
# the gamma distribution plot.
# What kind of parameters are these?

# Increasing alpha:
for (i in seq(5,25,5))
{
  x = seq(0, 10, length=200);
  y=dgamma(x,i,2) 
  plot(x, y, type="l", lwd=2, col="red");
}

# Increasing lamda = 1/beta:
i = 0;
for (i in seq(2,30,5))
{
  x = seq(0, 10, length=200);
  y=dgamma(x,5,i) 
  plot(x, y, type="l", lwd=2, col="red");
}

# A: It is obvious that changing alpha will result in a shift of the graph
# (and a slight scale). Changing lamda will result in a scale of the graph
# (and a slight shift).
# Alpha is the location parameter, and beta is the shape parameter.

# ******** TESTING QUANTILES ********:
qgamma(0.95, 5, 2)

## [1] 4.57676

#Hence, there is a 95% probability that a random number less than or equal to
# 4.577 is chosen from the gamma distribution. 

qgamma(0.20, 5,2)

## [1] 1.54477

# We now know that the probability of selecting a number from the 
# gamma distribution that is greater than or equal to 1.545 is 0.80. 

#:::::::::::::::::::::::END 1-GAMMA DISTRIBUTION:::::::::::::::::::::::::::::::::#

#:::::::::::::::::::::::1-EXPONENTIAL DISTRIBUTION::::::::::::::::::::::::::::::::::::#
# ******** THE PROBABILITY DENSITY FUNCTION ********:

# lamda = 2 
x = seq(0, 10, length=200);
y=dexp(x,2) 
plot(x, y, type="l", lwd=2, col="red");

# Q: Check the effect of changing the parameter and notice the effects on
# the exponential distribution plot.
# What kind of parameters are these?

# Increasing lamda = 1/beta:
i = 0;
for (i in seq(2,30,5))
{
  x = seq(0, 10, length=200);
  y=dexp(x,i) 
  plot(x, y, type="l", lwd=2, col="red");
}

# A: It is obvious that we cannot change alpha (no shift)
# But changing lamda will result in a scale of the graph (increase lambda, the graph becomes steeper)
# Thus lambda is the shape parameter.

# ******** TESTING QUANTILES ********:
qexp(0.95, 2)

## [1] 1.497866

#Hence, there is a 95% probability that a random number less than or equal to
# 1.498 is chosen from this distribution. 

qexp(0.20,2)

## [1] 0.1115718

# We now know that the probability of selecting a number from this
# distribution that is greater than or equal to .112 is 0.80. 

#:::::::::::::::::::::::END 1-EXPONENTIAL DISTRIBUTION:::::::::::::::::::::::::::::::::#

#:::::::::::::::::::::::1-CHI-SQUARED DISTRIBUTION::::::::::::::::::::::::::::::::::::#
# ******** THE PROBABILITY DENSITY FUNCTION ********:

# x= seq(0, 10, length=200)
# alpha = 5; 
# lamda = 2;
# y= (x^(alpha-1)*exp(-lamda*x)*(lamda)^alpha / gamma(alpha));
# plot(x, y, type="l", lwd=2, col="red") 
# 
# # is the same as:

x = seq(0, 10, length=200);
y=dchisq(x, df=1) 
plot(x, y, type="l", lwd=2, col="red");

# Q: Check the effect of changing the parameter and notice the effects on
# the chi-squared distribution plot.
# What kind of parameters are these?

# Increasing df (degrees of freedom):
for (i in seq(1,25,5))
{
  x = seq(0, 10, length=200);
  y=dchisq(x, df=i) 
  plot(x, y, type="l", lwd=2, col="red");
}

# A: Increasing the degrees of freedom will cause both a scale and a shift. It is a location
# and scale parameter. 

# ******** TESTING QUANTILES ********:
qchisq(0.95, 1)

## [1] 3.841459

#Hence, there is a 95% probability that a random number less than or equal to
# 3.841 is chosen from the standard normal distribution. 

qchisq(0.20, 1)

## [1] 0.06418475

# We now know that the probability of selecting a number from the standard 
# normal distribution that is greater than or equal to 0.0642 is 0.80. 

#:::::::::::::::::::::::END 1-CHI-SQUARED DISTRIBUTION:::::::::::::::::::::::::::::::::#

#:::::::::::::::::::::::1-BETA DISTRIBUTION::::::::::::::::::::::::::::::::::::#
# ******** THE PROBABILITY DENSITY FUNCTION ********:

x=seq(0,1,length=200)
y=dbeta(x,3,3)
plot(x,y,type="l",lwd=2,col="blue")

# Q: Check the effect of changing the two parameters and notice the effects on
# this distribution plot.
# What kind  of parameters are these?

# Increasing shape1:
i=0
for (i in seq(3,10,1))
{
  x=seq(0,1,length=200)
  y=dbeta(x,i,3)
  plot(x,y,type="l",lwd=2,col="blue")
}

# Increasing shape2:
i=0
for (i in seq(3,10,1))
{
  x=seq(0,1,length=200)
  y=dbeta(x,3,i)
  plot(x,y,type="l",lwd=2,col="blue")
}

# A: It is obvious that increasing shape1 causes a skewed-left graph and 
# increasing shape2 causes a skewed-right graph. They both shift and scale.
# Shape1 is the skewed-left location/shape parameter, and beta is the 
# skewed-right location/shape parameter.

# ******** TESTING QUANTILES ********:
qbeta(0.95,3,3)

## [1] 0.8107446

#Hence, there is a 95% probability that a random number less than or equal to
# .811 is chosen from this distribution. 

qbeta(0.20,3,3)

## [1] 0.3265979

# We now know that the probability of selecting a number from this 
#  distribution that is greater than or equal to .327 is 0.80. 

#:::::::::::::::::::::::END 2-BETA DISTRIBUTION:::::::::::::::::::::::::::::::::#

#:::::::::::::::::::::::2-UNIFORM DISTRIBUTION::::::::::::::::::::::::::::::::::::#
# ******** THE PROBABILITY DENSITY FUNCTION ********:

x=seq(0,1,length=200)
y=dunif(x, 0, 1)
plot(x,y,type="l",lwd=2,col="blue")

# Q: Check the effect of changing the two parameters and notice the effects on
# this distribution plot.
# What kind  of parameters are these?

# Increasing min:
i=0
for (i in seq(-25,0,5))
{
  x=seq(0,1,length=200)
  y=dunif(x, i, 1)
  plot(x,y,type="l",lwd=2,col="blue",ylim=c(0,3))
}

# Increasing max:
i=0
for (i in seq(1,25,5))
{
  x=seq(0,1,length=200)
  y=dunif(x, 0, i)
  plot(x,y,type="l",lwd=2,col="blue",ylim=c(0,3))
}

# A: It is obvious that increasing our min causes an enlarged scale, while  
# increasing our max causes a decreased scale. The graph remains between 0 and 1.
# Thus, the min and max are both scale parameters.

# ******** TESTING QUANTILES ********:
qunif(0.95,0,1)

## [1] 0.95

#Hence, there is a 95% probability that a random number less than or equal to
# 0.95 is chosen from this distribution. 

qbeta(0.20,0,1)

## [1] 0

# We now know that the probability of selecting a number from this 
#  distribution that is greater than or equal to 0 is 0.80. 

#:::::::::::::::::::::::END 2-UNIFORM DISTRIBUTION:::::::::::::::::::::::::::::::::#

#:::::::::::::::::::::::3-CAUCHY DISTRIBUTION::::::::::::::::::::::::::::::::::::#
# ******** THE PROBABILITY DENSITY FUNCTION ********:

x=seq(0,10,length=200)
y=dcauchy(x,100, 1)
plot(x,y,type="l",lwd=2,col="blue")

x=seq(0,10,length=200)
y=dcauchy(x,1/2, 1)
plot(x,y,type="l",lwd=2,col="blue")

# Q: Check the effect of changing the two parameters and notice the effects on
# this distribution plot.
# What kind  of parameters are these?

# Increasing first parameter:
i=0
for (i in seq(1/2,5,1))
{
  x=seq(0,10,length=200)
  y=dcauchy(x,i, 1)
  plot(x,y,type="l",lwd=2,col="blue",ylim=c(0,3))
}

# Increasing second parameter:
i=0
for (i in seq(1,5,1))
{
  x=seq(0,10,length=200)
  y=dcauchy(x,1/2, i)
  plot(x,y,type="l",lwd=2,col="blue",ylim=c(0,3))
}

# A: It is obvious that increasing our first parameter causes shift, while  
# increasing our second parameter causes a scale.
# Thus, the first parameter is a location parameter, while our second
# parameter is a scale parameter

# ******** TESTING QUANTILES ********:
qcauchy(0.95,100)

## [1] 106.3138

#Hence, there is a 95% probability that a random number less than or equal to
# 106.3138 is chosen from this distribution. 
qcauchy(0.95,1/2)

## [1] 6.813752

#Hence, there is a 95% probability that a random number less than or equal to
# 6.814 is chosen from this distribution. 


qcauchy(0.20,100)

## [1] 98.62362

# We now know that the probability of selecting a number from this 
# distribution that is greater than or equal to 98.62 is 0.80. 

qcauchy(0.20,1/2)

## [1] -0.8763819

# We now know that the probability of selecting a number from this 
# distribution that is greater than or equal to 0.8763 is 0.80. 
#:::::::::::::::::::::::END 3-CAUCHY DISTRIBUTION:::::::::::::::::::::::::::::::::#

#:::::::::::::::::::::::4-WEIBULL DISTRIBUTION::::::::::::::::::::::::::::::::::::#
# ******** THE PROBABILITY DENSITY FUNCTION ********:

x=seq(0,10,length=200)
y=dweibull(x,3,3)
plot(x,y,type="l",lwd=2,col="blue")

# Q: Check the effect of changing the two parameters and notice the effects on
# this distribution plot.
# What kind  of parameters are these?

# Increasing first parameter:
i=0
for (i in seq(3,15,3))
{
  x=seq(0,10,length=200)
  y=dweibull(x,i,3)
  plot(x,y,type="l",lwd=2,col="blue",ylim=c(0,3))
}

# Increasing second parameter:
i=0
for (i in seq(3,10,1.5))
{
  x=seq(0,10,length=200)
  y=dweibull(x,3,i)
  plot(x,y,type="l",lwd=2,col="blue",ylim=c(0,3))
}

# A: It is obvious that increasing our first parameter causes a change in shape, while  
# increasing our second parameter causes a shift/scale
# Thus, the first parameter is a shape parameter, while our second
# parameter is a scale parameter

# ******** TESTING QUANTILES ********:
qweibull(0.95,3,3)

## [1] 4.324696

#Hence, there is a 95% probability that a random number less than or equal to
# 4.32 is chosen from this distribution. 

qweibull(0.20,3,3)

## [1] 1.819628

# We now know that the probability of selecting a number from this 
# distribution that is greater than or equal to 1.820 is 0.80. 

#:::::::::::::::::::::::END 4-WEIBULL DISTRIBUTION:::::::::::::::::::::::::::::::::#

# PART (B) - CHECKING EMPIRICAL RULES ON EACH DISTRIBUTION:

# GAMMA DISTRIBUTION
# ******** 68%-95%-99.7% Rule {Empirical Rule} ********:
# 1 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rgamma(10000,5,2))
standd = sqrt(var(rgamma(10000,5,2)))
numDeviations = 1
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dgamma(x,5,2) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dgamma(x,5,2) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pgamma(mu + standdScale,5,2) - pgamma(mu - standdScale,5,2) 

## [1] 0.6935619

# 2 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rgamma(10000,5,2))
standd = sqrt(var(rgamma(10000,5,2)))
numDeviations = 2
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dgamma(x,5,2) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dgamma(x,5,2) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pgamma(mu + standdScale,5,2) - pgamma(mu - standdScale,5,2) 

## [1] 0.9580092

# 3 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rgamma(10000,5,2))
standd = sqrt(var(rgamma(10000,5,2)))
numDeviations = 3
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dgamma(x,5,2) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dgamma(x,5,2) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pgamma(mu + standdScale,5,2) - pgamma(mu - standdScale,5,2) 

## [1] 0.9909296

# Looking at the above for this distribution, we see that it almost worked, but
# not just quite. So no.
# 1sd = 70.8%
# 2sd = 95.5%
# 3sd = 99.05%

# EXPONENTIAL DISTRIBUTION
# ******** 68%-95%-99.7% Rule {Empirical Rule} ********:
# 1 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rexp(10000,2))
standd = sqrt(var(rexp(10000,2)))
numDeviations = 1
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dexp(x,2) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dexp(x,2) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pexp(mu + standdScale,2) - pexp(mu - standdScale,2) 

## [1] 0.8648885

# 2 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rexp(10000,2))
standd = sqrt(var(rexp(10000,2)))
numDeviations = 2
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dexp(x,2) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dexp(x,2) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pexp(mu + standdScale,2) - pexp(mu - standdScale,2) 

## [1] 0.9489496

# 3 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rexp(10000,2))
standd = sqrt(var(rexp(10000,2)))
numDeviations = 3
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dexp(x,2) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dexp(x,2) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pexp(mu + standdScale,2) - pexp(mu - standdScale,2) 

## [1] 0.980742

# Looking at the above for this distribution, we see that NO it does not work.
# 1sd = 85%
# 2sd = 95.0%
# 3sd = 97.99%

# CHI SQUARED DISTRIBUTION
# ******** 68%-95%-99.7% Rule {Empirical Rule} ********:
# 1 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rchisq(10000,df=1))
standd = sqrt(var(rchisq(10000,df=1)))
numDeviations = 1
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dchisq(x,df=1) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dchisq(x,df=1) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pchisq(mu + standdScale,1) - pchisq(mu - standdScale,1) 

## [1] 0.879188

# 2 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rchisq(10000,df=1))
standd = sqrt(var(rchisq(10000,df=1)))
numDeviations = 2
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dchisq(x,df=1) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dchisq(x,df=1) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pchisq(mu + standdScale,1) - pchisq(mu - standdScale,1) 

## [1] 0.9483769

# 3 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rchisq(10000,df=1))
standd = sqrt(var(rchisq(10000,df=1)))
numDeviations = 3
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dchisq(x,df=1) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dchisq(x,df=1) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pchisq(mu + standdScale,1) - pchisq(mu - standdScale,1) 

## [1] 0.9782385

# Looking at the above for this distribution, we see that NO it does not work.
# 1sd = 88.1%
# 2sd = 94.98%
# 3sd = 97.83%

# BETA DISTRIBUTION
# ******** 68%-95%-99.7% Rule {Empirical Rule} ********:
# 1 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rbeta(10000,3,3))
standd = sqrt(var(rbeta(10000,3,3)))
numDeviations = 1
standdScale = numDeviations * standd

x = seq(0, 1, length=200);
y=dbeta(x,3,3) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dbeta(x,3,3) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pbeta(mu + standdScale,3,3) - pbeta(mu - standdScale,3,3) 

## [1] 0.6388164

# 2 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rbeta(10000,3,3))
standd = sqrt(var(rbeta(10000,3,3)))
numDeviations = 2
standdScale = numDeviations * standd

x = seq(0, 1, length=200);
y=dbeta(x,3,3) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dbeta(x,3,3) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pbeta(mu + standdScale,3,3) - pbeta(mu - standdScale,3,3) 

## [1] 0.9735811

# 3 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rbeta(10000,3,3))
standd = sqrt(var(rbeta(10000,3,3)))
numDeviations = 3
standdScale = numDeviations * standd

x = seq(0, 1, length=200);
y=dbeta(x,3,3) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dbeta(x,3,3) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pbeta(mu + standdScale,3,3) - pbeta(mu - standdScale,3,3) 

## [1] 1

# Looking at the above for this distribution, we see that NO it does not work.
# 1sd = 64.6%
# 2sd = 96.7%
# 3sd = 100%

# UNIFORM DISTRIBUTION
# ******** 68%-95%-99.7% Rule {Empirical Rule} ********:
# 1 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(runif(10000,0,1))
standd = sqrt(var(runif(10000,0,1)))
numDeviations = 1
standdScale = numDeviations * standd

x = seq(0, 1, length=200);
y=dunif(x,0,1) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dunif(x,0,1) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

punif(mu + standdScale,0,1) - punif(mu - standdScale,0,1) 

## [1] 0.5805254

# 2 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(runif(10000,0,1))
standd = sqrt(var(runif(10000,0,1)))
numDeviations = 2
standdScale = numDeviations * standd

x = seq(0, 1, length=200);
y=dunif(x,0,1) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dunif(x,0,1) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

punif(mu + standdScale,0,1) - punif(mu - standdScale,0,1) 

## [1] 1

# 3 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(runif(10000,0,1))
standd = sqrt(var(runif(10000,0,1)))
numDeviations = 3
standdScale = numDeviations * standd

x = seq(0, 1, length=200);
y=dunif(x,0,1) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dunif(x,0,1) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

punif(mu + standdScale,0,1) - punif(mu - standdScale,0,1) 

## [1] 1

# Looking at the above for this distribution, we see that NO it does not work.
# 1sd = 58.4%
# 2sd = 100%
# 3sd = 100%

# CAUCHY DISTRIBUTION
# ******** 68%-95%-99.7% Rule {Empirical Rule} ********:
# 1 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rcauchy(10000,100,1))
standd = sqrt(var(rcauchy(10000,100,1)))
numDeviations = 1
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dcauchy(x,100,1) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dcauchy(x,100,1) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pcauchy(mu + standdScale,100,1) - pcauchy(mu - standdScale,100,1) 

## [1] 0.9952891

# 2 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rcauchy(10000,100,1))
standd = sqrt(var(rcauchy(10000,100,1)))
numDeviations = 2
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dcauchy(x,100,1) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dcauchy(x,100,1) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pcauchy(mu + standdScale,100,1) - pcauchy(mu - standdScale,100,1) 

## [1] 0.9960959

# 3 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rcauchy(10000,100,1))
standd = sqrt(var(rcauchy(10000,100,1)))
numDeviations = 3
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dcauchy(x,100,1) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dcauchy(x,100,1) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pcauchy(mu + standdScale,100,1) - pcauchy(mu - standdScale,100,1) 

## [1] 0.995258

# Looking at the above for this distribution, we see that NO it does not work.
# 1sd = 98.21%
# 2sd = 99.81%
# 3sd = 99.99%

# WEIBULL DISTRIBUTION
# ******** 68%-95%-99.7% Rule {Empirical Rule} ********:
# 1 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rweibull(10000,3,3))
standd = sqrt(var(rweibull(10000,3,3)))
numDeviations = 1
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dweibull(x,3,3) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dweibull(x,3,3) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pweibull(mu + standdScale,3,3) - pweibull(mu - standdScale,3,3) 

## [1] 0.6674745

# 2 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rweibull(10000,3,3))
standd = sqrt(var(rweibull(10000,3,3)))
numDeviations = 2
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dweibull(x,3,3) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dweibull(x,3,3) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pweibull(mu + standdScale,3,3) - pweibull(mu - standdScale,3,3) 

## [1] 0.9625031

# 3 SD From mean:
# First find mean, standard deviation, and set scale
mu = mean(rweibull(10000,3,3))
standd = sqrt(var(rweibull(10000,3,3)))
numDeviations = 3
standdScale = numDeviations * standd

x = seq(0, 10, length=200);
y=dweibull(x,3,3) 
plot(x, y, type="l", lwd=2, col="blue");
x=seq(mu - standdScale,mu + standdScale,length=200)
y=dweibull(x,3,3) 
polygon(c(mu - standdScale,x,mu + standdScale),c(0,y,0),col="gray") 

pweibull(mu + standdScale,3,3) - pweibull(mu - standdScale,3,3)  

## [1] 0.9986898

# Looking at the above for this distribution, we see that althought it is close,
# it does not work.
# 1sd = 66.72%
# 2sd = 96.17%
# 3sd = 99.86%


#PART C: Take a random sample of size 2000 from each distribution.
# Calculate its mean and standard deviation. Does is match the closed 
# form of Expected value and Variance derived in class

# GAMMA
sample = rgamma(2000, 5, 2)
summary(sample)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.2725  1.6680  2.3770  2.5050  3.1960  7.0340

var(sample)

## [1] 1.269144

sqrt(var(sample))

## [1] 1.126563

mean=5 * .5; mean

## [1] 2.5

var=5 * .5^2; var

## [1] 1.25

stddev=sqrt(var); stddev

## [1] 1.118034

# Yes. We see that the mean and variance closely match the expected value
# and variance equations we derived in class.

# EXPONENTIAL
sample = rexp(2000, 2)
summary(sample)

##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 0.000044 0.137400 0.330700 0.485200 0.656300 3.461000

var(sample)

## [1] 0.2501652

sqrt(var(sample))

## [1] 0.5001652

mean=1 / 2; mean

## [1] 0.5

var=1 / (2^2); var

## [1] 0.25

stddev=sqrt(var); stddev

## [1] 0.5

# Yes. We see that the mean and variance closely match the expected value
# and variance equations we derived in class.

# CHI-SQUARED
sample = rchisq(2000, df=1)
summary(sample)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.1112  0.4473  1.0050  1.3330 13.2200

var(sample)

## [1] 2.062137

sqrt(var(sample)) #standard dev

## [1] 1.436014

mean=1; mean #b/c E[x] = n = 1

## [1] 1

var=2; var #b/c Var[x] = 2n = 2

## [1] 2

stddev=sqrt(var); stddev

## [1] 1.414214

# Yes. We see that the mean and variance closely match the expected value
# and variance equations we derived in class.

# BETA
sample = rbeta(2000, 3,3)
summary(sample)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.03385 0.36380 0.50360 0.50410 0.64500 0.96170

var(sample)

## [1] 0.03590061

sqrt(var(sample)) #standard dev

## [1] 0.1894746

mean=3/(3+3); mean #b/c E[x] = a/(a+b)

## [1] 0.5

var=(3*3)/((3+3)^2*(3+3+1)); var #b/c Var[x] = ab/((a+b)^2 (a+b+1))

## [1] 0.03571429

stddev=sqrt(var); stddev

## [1] 0.1889822

# Yes. We see that the mean and variance closely match the expected value
# and variance equations we derived in class.

# UNIFORM
sample = runif(2000, 0,1)
summary(sample)

##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 0.0000342 0.2389000 0.4947000 0.4942000 0.7487000 0.9999000

var(sample)

## [1] 0.08493495

sqrt(var(sample)) #standard dev

## [1] 0.291436

mean=(1+0)/2; mean #b/c E[x] = (a+b)/2

## [1] 0.5

var=(1-0)^2/12; var #b/c Var[x] = (b-a)^2/12

## [1] 0.08333333

stddev=sqrt(var); stddev

## [1] 0.2886751

# Yes. We see that the mean and variance closely match the expected value
# and variance equations we derived in class.

# # CAUCHY
# sample = rcauchy(2000, 100,1)
# summary(sample)
# var(sample)
# sqrt(var(sample)) #standard dev
# 
# mean=
# var=
# stddev=sqrt(var); stddev
# 
# sample = rcauchy(2000, 1/2,1)
# summary(sample)
# var(sample)
# sqrt(var(sample)) #standard dev
# 
# mean=
# var=
# stddev=sqrt(var); stddev
# 
# 
# # ??? We see that the mean and variance closely match the expected value
# # and variance equations we derived in class (IN BOTH CASES).
# 
# # WEIBULL
# sample = rweibull(2000, 3,3)
# summary(sample)
# var(sample)
# sqrt(var(sample)) #standard dev
# 
# mean=gamma(1/3 + 1); mean #b/c E[x] = Gamma(1/alpha + 1)
# var=gamma(2/3 + 1) - (gamma(1/3 + 1))^2; var #b/c Var[x] = Gamma(2/alpha + 1) - Gamma^2(1/alpha +1)
# stddev=sqrt(var); stddev
# 
# # Yes. We see that the mean and variance closely match the expected value
# # and variance equations we derived in class.